Examen parcial 2 Estadística No Paramétrica
Marcelino Sánchez Rodríguez
191654
Pregunta 1
Usando los datos FishMercury del paquete resampledata, que contiene los niveles de mercurio (en partes por millón, ppm) para 30 peces capturados en lagos de Minnesota: Hacer un histograma de los datos. Describir lo que se observa. Hacer 200 réplicas bootstrap de la media y obtener el error estándar y el intervalo bootstrap percentil 95 %.
## cargamos paquetes necesarios y data
library(resampledata)##
## Attaching package: 'resampledata'## The following object is masked from 'package:datasets':
##
## Titanicdata1 <- FishMercury$Mercury
hist(FishMercury[,1], main = "Histograma de los datos", xlab = "Niveles de mercurio (ppm)", ylab = "Frecuencia", col = "blue", breaks= 25)
Observamos que la mayoría de los datos están concentrados cerca del 0 y al parecer hay un outlier cerca de 2 ppm.
Quitar el outlier, y hacer réplicas boostrap de la media de los datos sin el outlier. Obtener el error estándar y el intervalo bootstrap percentil 95 %.
library(boot)
library(bootstrap)
set.seed(100)
## sin quitar outlier
Tn <- function(data, index){mean(data[index])}
ndata <- length(data1)
bootsFish2 <- boot(data = data1, statistic = Tn, R = 200)
(mean(bootsFish2$t) ) # estimación de media## [1] 0.1842095(sd(bootsFish2$t)) #estimador de desviación estándar## [1] 0.05828282#intervalo de confianza
(confFish <- boot.ci(bootsFish2, conf = .95, type = "perc"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 200 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bootsFish2, conf = 0.95, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 0.1098, 0.3067 )
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable#distribución bootstrap
hist(bootsFish2$t, breaks = 50)
##---------------------------------------------------------------------------------------------
##quitando outlier
bootsFish3 <- boot(data = data1[-1], statistic = Tn, R = 200)
(mean(bootsFish3$t) ) # estimación de media## [1] 0.1239441(sd(bootsFish3$t)) #estimador de desviación estándar## [1] 0.00845829#intervalo de confianza
(confFish2 <- boot.ci(bootsFish3, conf = .95, type = "perc"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 200 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bootsFish3, conf = 0.95, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 0.1048, 0.1401 )
## Calculations and Intervals on Original Scale
## Some percentile intervals may be unstable#distribución bootstrap
hist(bootsFish3$t, breaks = 50)
¿Qué efecto tuvo el remover el outlier en la distribución bootstrap, y en el error estándar?
Hizo que el intervalo de confianza fuera más chica ya que la desviación estándar fue menor, y que la distribución bootstrap fuera más simétrica.
Pregunta 2
- Con los datos Girls2004 del paquete resampledata,
library(resampledata)
data2 <- Girls2004Obtener estadísticas resumen de los pesos de las bebés nacidas en Wyoming (WY) y Alaska (AK) (hacer un análisis separado para cada estado).
# Estadísticas resumen para bebés nacidos en Wyoming (WY)
wy_data <- subset(data2, State == "WY")
summary(wy_data$Weight)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2212 2934 3278 3208 3515 3995hist(wy_data$Weight, main = "Histograma de los datos WY", xlab = "Peso (libras)", ylab = "Frecuencia", col = "blue", breaks= 15)
# Estadísticas resumen para bebés nacidos en Alaska (AK)
ak_data <- subset(data2, State == "AK")
summary(ak_data$Weight)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2182 3170 3558 3516 3926 4592hist(ak_data$Weight, main = "Histograma de los datos AK", xlab = "Peso (libras)", ylab = "Frecuencia", col = "blue", breaks= 15)
Hacer bootstrap de la diferencia de medias, graficar la distribución y dar las estadísticas sumarias. Obtener un intervalo boostrap percentil de 95 %, e interpretar este intervalo.
Bnum <- 1000
wy_list <- wy_data$Weight
ak_list <- ak_data$Weight
diff_boot_list <- c()
for (i in 1:Bnum) {
wy_sample <- sample(wy_list, length(wy_list), replace = TRUE)
ak_sample <- sample(ak_list, length(ak_list), replace = TRUE)
wy_boot <- mean(wy_sample)
ak_boot <- mean(ak_sample)
diff_boot <- wy_boot - ak_boot
diff_boot_list <- c(diff_boot_list, diff_boot)
}
hist(diff_boot_list, main = "Histograma de la diferencia de medias", xlab = "Diferencia de medias", ylab = "Frecuencia", col = "blue", breaks= 15)
summary(diff_boot_list)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -725.9 -387.5 -309.6 -310.1 -235.6 170.2#intervalo de confianza
diff_boot_listOrdered <- sort(diff_boot_list)
lower <- diff_boot_listOrdered[round(0.025 * Bnum)]
upper <- diff_boot_listOrdered[round(0.975 * Bnum)]
lower## [1] -535.65upper## [1] -93.175En primer lugar dada la distribución notamos que en general los pesos de los bebés nacidos en Wyoming son menores que los de los bebés nacidos en Alaska. El intervalo de confianza nos dice que la diferencia de medias se encuentra entre 0.5 y 0.1 libras, lo que significa que en promedio los bebés nacidos en Wyoming pesan 0.3 libras menos que los bebés nacidos en Alaska.
Obtener el estimador del sesgo bootstrap. ¿Qué fracción representa del error estándar bootstrap?
estDiffReal <- mean(wy_list) - mean(ak_list)
sesgo <- - mean(diff_boot_list)-estDiffReal
sesgo/(sd(diff_boot_list))# checar qué onda con esto## [1] 5.430541Realizar una prueba de permutación para ver si la diferencia de medias en los pesos es significativa estadísticamente y dar una conclusión.
Bnum <- 1000
z <- c(wy_list, ak_list)
ind<- 1:length(z)
reps <- numeric(Bnum)
t0 <- t.test(wy_list, ak_list)$statistic
for (i in 1:Bnum) {
indices <- sample(ind, size = length(wy_list), replace = F)
x1 <- z[indices]
x2 <- z[-indices]
reps[i] <- t.test(x1,x2)$statistic
}
p <- mean(c(t0,reps) > t0)
p## [1] 0.995005No es significativa, por lo que no podemos rechazar la hipótesis nula de que las medias de los pesos de los bebés nacidos en Wyoming y Alaska son iguales. Es decir, en realidad no hay diferencias en tre los pesos de los distintos estados.
# Pregunta 3
La función smooth.spline del paquete básico stats permite dar los grados de libertad explícitamente. Escribir una función que elija el número de grados de libertad vía validación cruzada.
```r
library(stats)
cv_smooth_spline <- function(x, y, nk = 10, df_seq = seq(2, 20, by = 2)) {
# nk: número de particiones para la validación cruzada
# df_seq: secuencia de grados de libertad a probar
# datos
n <- length(x)
# indices de los datos
idx <- seq_len(n)
# particiones para la validación cruzada
folds <- sample(rep(1:nk, length.out = n))
# mse para cada df
mse <- rep(NA, length(df_seq))
for (i in seq_along(df_seq)) {
# df actual
df <- df_seq[i]
# mse para cada partición
cv_mse <- rep(NA, nk)
for (k in seq_len(nk)) {
# indices de los datos de entrenamiento
train_idx <- which(folds != k)
# ajustar el modelo
fit <- smooth.spline(x[idx[train_idx]], y[idx[train_idx]], df = df)
# mse en la partición k
cv_mse[k] <- mean((y[idx[-train_idx]] - predict(fit, x[idx[-train_idx]])$y)^2)
}
# mse promedio para el df actual
mse[i] <- mean(cv_mse)
}
# df con el mse más pequeño
best_df <- df_seq[which.min(mse)]
# regresar los resultados
return(list(best_df = best_df, mse = mse, df_seq = df_seq))
}
## EJEMPLO
set.seed(191654)
x <- c(runif(100, 0, 10))
y <- c(cos(x) + rnorm(100, 0, 0.5))
k <- 5
plot(x, y)
# Estimar el número óptimo de grados de libertad usando validación cruzada
salida <- cv_smooth_spline(x, y, nk = k, df_seq = seq(2, 20, by = 1))
salida[1]## $best_df
## [1] 7# Ajustar el smooth spline con el número óptimo de grados de libertad
fit <- smooth.spline(x, y, df=salida[1])
# Graficar los datos y la curva ajustada
plot(x, y, pch=20, col="blue", xlab="x", ylab="y")
lines(fit, col="red")
legend("topright", legend=c("datos", "smooth spline"), col=c("blue", "red"), pch=20, lty=1)
Pregunta 4
El archivo earthquake.dat contiene registros de niveles de agua para un conjunto de seis pozos en California (descripción detallada: earthquake.txt.Las medidas se hicieron a lo largo del tiempo. Cons- truir un estimador LOESS para examinar tendencias en los datos. ¿El estimador LOESS es exitoso? ¿En dónde falla para capturar las tendencias de los datos?
# Cargamos datos
data4 <- read.fwf("https://raw.githubusercontent.com/jvega68/ENP/main/datos/earthquake.dat",
width = c(17,7,1,7,1,6,4,2,2,80),
header = F, sep ="\t")
names(data4) <- c("fecha","lat","lat2", "lon","lon2", "depth","mag","lt","Q", "lugar")
t <- 1:length(data4$fecha)
plot(t,data4$mag)
# Estimador LOESS
loessFit <- loess(data4$mag ~ t, span = 0.5)
grid_span <- c(0.2,0.3,0.4,0.5,0.6,0.7)
# Graficar varios bandwidths
for(h in grid_span){
m1 <- loess(data4$mag ~ t, span = h, degree = 2)
lines(fitted(m1), col = which(grid_span == h), lwd = 3)
}
legend("topright",col = 1:6, legend = grid_span[1:6], lty = 1, lwd = 3)
Por la posición de los datos está ajustando una línea horizontal en todos los distintos valores del span que definimos para LOESS. Y en cierta forma, me parece qu está siendo un buen estimador, dado que todos los datos están posicionados de forma horizontal. Sin embargo, igual se nota que en los extremos el estimador LOESS deja de comportarse como una línea horizontal y se dispersa a disitintas direcciones según el valor del span. Repetir el ejercicio con NW y con polinomios locales.
# Nadarayan-Watson
library(locfit)## locfit 1.5-9.7 2023-01-02alpha <- seq(0.15, 0.99, by = 0.01)
a <- gcvplot(data4$mag ~ t, kern = "gauss", deg = 0, data = data4, alpha = alpha)
plot(a)
# Valor de alpha que minimiza los valores de a
alfa_star <- a$alpha[which.min(a$values)]
alfa_star## [1] 0.15# Modelo con alpha_star
modelo <- locfit(data4$mag ~ t, kern = "gauss", deg = 0, data = data4, alpha = alfa_star)
summary(modelo)## Estimation type: Local Regression
##
## Call:
## locfit(formula = data4$mag ~ t, data = data4, kern = "gauss",
## deg = 0, alpha = alfa_star)
##
## Number of data points: 2476
## Independent variables: t
## Evaluation structure: Rectangular Tree
## Number of evaluation points: 31
## Degree of fit: 0
## Fitted Degrees of Freedom: 9.356#comparando distintos modelos con distinta alpha
b0 <- locfit(data4$mag ~ t, kern="gauss", data = data4, deg = 0, alpha = 0.15)
b9 <- locfit(data4$mag ~ t, kern="gauss", data = data4, deg = 0, alpha = 0.99)
m2 <- locfit(data4$mag ~ t, kern="gauss", data = data4, deg = 0, alpha = alfa_star)
# Gráfica de los modelos
plot(data4$mag , col="darkgrey", main = "Nadaraya-Watson")
legend("bottomleft",
lty = 1,
col = c("red","blue","orange"),
legend=c(paste("bandwidth =", c( 0.15, 0.99, alfa_star))), bty = "n", lwd = 1.5)
lines(fitted(b0), col = "red", lwd = 1.5)
lines(fitted(b9), col = "blue", lwd = 1.5)
lines(fitted(m2), col = "orange", lwd = 3) 
# calcula ic de 95% con estimador local de varianza
crit(m2) <- crit(m2, cov = 0.95)
plot(m2, band = "local") # gráfica
predict(m2, se = T, band= "local") # valores## $fit
## [1] 4.715388 4.504636 4.395476 4.808349 4.711048 4.668345 4.684775 4.734160
## [9] 4.712907 4.811893 4.708274 4.708724 4.729118 4.706921 4.621296 4.711893
## [17] 4.466512 4.514456 4.313763 4.496028 4.468812 4.313644 4.520156 4.471362
## [25] 4.526623 4.571347 4.507701 4.523367 4.502347 4.600009 4.630093
##
## $se.fit
## [1] 0.05750472 0.05302356 0.06314133 0.05558678 0.05331504 0.05446211
## [7] 0.05552311 0.05712480 0.05395619 0.05677660 0.05158681 0.05117624
## [13] 0.05496119 0.04985991 0.05662270 0.05381084 0.05914577 0.05818672
## [19] 0.06740254 0.06179859 0.06371979 0.06588064 0.05472861 0.06175891
## [25] 0.05388179 0.04823814 0.05599660 0.05990951 0.05105558 0.04800379
## [31] 0.04804102
##
## $residual.scale
## rv
## 0.912385# Ajuste de polinomios locales de grado 3 y 4
pol_local_3 <- locfit(data4$mag ~ t, deg = 3, data = data4)
plot(data4$mag, lwd=2)
lines(pol_local_3, col = "red", lwd = 2)
#--------------------------------------------
pol_local_4 <- locfit(data4$mag ~ t, deg = 4, data = data4)
lines(pol_local_4, col = "blue", lwd = 2)
En estos nuevos estimadores notamos que se sigue una tendencia de lineal horizontal. Sin embargo, en estos casos notamos que en el estimador NW hay una mejor estabilidad en los extremos y de cierta forma logra encontrar cierta dinámica que ocurren en el centro de la gráfica donde hay mayor cantidad de datos. # Problema 5
Simular un conjunto de datos como sigue:
set.seed(191654)
x <- sort(runif(100))
y <- x^2 + 0.1*rnorm(100)Ajustar un spline a los datos simulados.
# Ajustar un spline a los datos simulados
splineFit <- smooth.spline(x,y, df = 5)
plot(x,y, pch = 19, col = "blue")
lines(x, predict(splineFit)$y, col = "red", lwd = 2)
# Ajustamos un spline de interpolación a los datos simulados
fit <- splinefun(x,y)
xx <- seq(0, max(x), length = 100)
yy <- fit(xx)
plot(x, y)
lines(xx, yy, lwd=3, col="green")
Problema 6
Con los datos en el archivo motor.dat, utilizando las variables time (tiempo en milisegundos) y accel (aceleración de la cabeza medida en g), ajustar los siguientes modelos: NW Polinomio local de grado 3 Spline cúbico Obtener en cada caso, la matriz sombrero. ¿Qué modelo es mejor?
library(locfit)
data6 <- read.table("/Users/marcelinosanchezrodriguez/itamRepositories/estNoPar2023Prim/examenEstNoPar/motor.dat", header = T, sep ="\t")# Ajuste del polinomio local de grado 3
pol_local_3 <- locfit(data6$accel ~ lp(data6$times, scale = T, nn = 0.5, deg = 3), data = data6)
plot(lp(data6$times,data6$accel))
lines(pol_local_3)
# Matriz Hat para el modelo de polinomio local de grado 3
#hat_pol_3 <- predict(diag(4), se.fit = TRUE, type = "hat")
#--------------------------------------------
# Ajuste del spline cúbico
splineFit <- smooth.spline(data6$times,data6$accel, df = 10)
plot(data6$times,data6$accel, pch = 19, col = "blue", main = "Spline cúbico")
lines(data6$times, predict(splineFit)$y, col = "red", lwd = 2)
# Matriz Hat para el modelo de spline cúbico
#hat_spline <- predict(diag(4), se.fit = TRUE, type = "hat")
#--------------------------------------------
# Ajuste del modelo NW
alpha <- seq(0.15, 0.99, by = 0.01)
a <- gcvplot(data6$accel ~ data6$times, kern = "gauss", deg = 0, data = data6, alpha = alpha)
# Valor de alpha que minimiza los valores de a
alfa_star <- a$alpha[which.min(a$values)]
alfa_star## [1] 0.19modelo <- locfit(data6$accel ~ data6$times, kern = "gauss", deg = 0, data = data6, alpha = alfa_star)
summary(modelo)## Estimation type: Local Regression
##
## Call:
## locfit(formula = data6$accel ~ data6$times, data = data6, kern = "gauss",
## deg = 0, alpha = alfa_star)
##
## Number of data points: 94
## Independent variables: data6$times
## Evaluation structure: Rectangular Tree
## Number of evaluation points: 23
## Degree of fit: 0
## Fitted Degrees of Freedom: 7.941b0 <- locfit(data6$accel ~ data6$times, kern="gauss", data = data6, deg = 0, alpha = 0.15)
b9 <- locfit(data6$accel ~ data6$times, kern="gauss", data = data6, deg = 0, alpha = 0.99)
m2 <- locfit(data6$accel ~ data6$times, kern="gauss", data = data6, deg = 0, alpha = alfa_star)
plot(data6$accel , col="darkgrey", main = "Nadaraya-Watson")
legend("bottomleft",
lty = 1,
col = c("red","blue","orange"),
legend=c(paste("bandwidth =", c( 0.15, 0.99, alfa_star))), bty = "n", lwd = 1.5)
lines(fitted(b0), col = "red", lwd = 1.5)
lines(fitted(b9), col = "blue", lwd = 1.5)
lines(fitted(m2), col = "orange", lwd = 3)
# Matriz Hat para el modelo NW
#hat_nw <- predict(diag(4), se.fit = TRUE, type = "hat")Problema 7
En el problema anterior, para cada uno de los modelos, Calcular una banda de confianza del 95 % alrededor del estimador.
# Banda de confianza del 95% para el modelo de polinomio local de grado 3
# Paso 1 -- gráfica de replicas individuales
plot(data6$times,data6$accel)
# Paso 2 -- Iteraciones bootstrap
set.seed(100)
alfa_star## [1] 0.19B <- 1000
yfit <- matrix(nrow = length(data6$times), ncol = B)
for (i in 1:B) {
ind <- sample(1:length(data6$times), size = length(data6$times), replace=TRUE)
# print(head(randsitenum))
x <- as.vector(data6$times[ind])
y <- as.vector(data6$accel[ind])
locboot <- locfit(y ~ lp(x, scale = T, nn = 0.5, deg = 3))
predboot <- predict(locboot, newdata = sort(x))
yfit[,i] <- predboot
# note plotting lines is slowww
lines(sort(x), yfit[,i], lwd=2, col = rgb(0.5,0.5,0.5,0.10))
}
# Paso 3 -- Datos ajustados originales
m4p <- predict(m2, newdata = data6$times)
lines(data6$times, m4p, lwd=2, col="blue")
# Step 4 -- Dibuja el intervalo bootstrap
yfit95 <- apply(yfit, 1, function(x) quantile(x, prob = 0.95, na.rm=T))
yfit05 <- apply(yfit, 1, function(x) quantile(x, prob = 0.05, na.rm=T))
lines(data6$times, yfit95, lwd=1, col="red")
lines(data6$times, yfit05, lwd=1, col="red")
#--------------------------------------------
# Banda de confianza del 95% para el modelo de spline cúbico
library(npreg)##
## Attaching package: 'npreg'## The following object is masked from 'package:boot':
##
## bootsplineFit <- ss(data6$times,data6$accel, df = 10)
a <- boot(splineFit, R = 10)##
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|============================ | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%plot(a)
#--------------------------------------------
# Banda de confianza del 95% para el modelo NW
# Paso 1 -- gráfica de replicas individuales
plot(data6$times,data6$accel)
# Paso 2 -- Iteraciones bootstrap
set.seed(100)
alfa_star## [1] 0.19B <- 1000
yfit <- matrix(nrow = length(data6$times), ncol = B)
for (i in 1:B) {
ind <- sample(1:length(data6$times), size = length(data6$times), replace=TRUE)
# print(head(randsitenum))
x <- as.vector(data6$times[ind])
y <- as.vector(data6$accel[ind])
locboot <- locfit(y ~ x, kern="gauss", deg = 0, alpha = alfa_star)
predboot <- predict(locboot, newdata = sort(x))
yfit[,i] <- predboot
# note plotting lines is slowww
lines(sort(x), yfit[,i], lwd=2, col = rgb(0.5,0.5,0.5,0.10))
}
# Paso 3 -- Datos ajustados originales
m4p <- predict(m2, newdata = data6$times)
lines(data6$times, m4p, lwd=2, col="blue")
# Step 4 -- Dibuja el intervalo bootstrap
yfit95 <- apply(yfit, 1, function(x) quantile(x, prob = 0.95, na.rm=T))
yfit05 <- apply(yfit, 1, function(x) quantile(x, prob = 0.05, na.rm=T))
lines(data6$times, yfit95, lwd=1, col="red")
lines(data6$times, yfit05, lwd=1, col="red")
Analizar si los residuales cumplen con los supuestos del modelo, verificando normalidad y varianza constante.
# Residuales del modelo de polinomio local de grado 3
resid_pol_3 <- resid(pol_local_3)
plot(resid_pol_3)
qqnorm(resid_pol_3)
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(resid_pol_3)##
## Shapiro-Wilk normality test
##
## data: resid_pol_3
## W = 0.98108, p-value = 0.1916#--------------------------------------------
# Residuales del modelo de spline cúbico
resid_spline <- resid(splineFit)
plot(resid_spline)
qqnorm(resid_spline)
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(resid_spline)##
## Shapiro-Wilk normality test
##
## data: resid_spline
## W = 0.97933, p-value = 0.1429#--------------------------------------------
# Residuales del modelo NW
resid_nw <- resid(m2)
plot(resid_nw)
qqnorm(resid_nw)
# Prueba de normalidad de Shapiro-Wilk
shapiro.test(resid_nw)##
## Shapiro-Wilk normality test
##
## data: resid_nw
## W = 0.98255, p-value = 0.2438En los 3 modelos notamos que tenemos problemas de homocedasticidad, por lo que tenemos fallos en los supuestos del modelo, en particular, se presenta en medio de la gráfica.
Por último, se puede notar por gráfica y test de Shapiro que los residuales siguen una distribución normal. En este caso se cumplen los supuestos.
Problema 8
Utilizando el conjunto de datos EuStockMarkets que es de los archivos de datos básicos de R, se tienen los precios diarios de 4 indicadores financieros: el índice DAX, el SMI, el CAC y el FTSE.
data8 <- EuStockMarketsCalcular los rendimientos diarios en logaritmos de los índices DAX y FTSE
# Seleccionar las columnas correspondientes a DAX y FTSE
dax_ftse <- data8[, c("DAX", "FTSE")]
# Calcular los rendimientos diarios en logaritmos de DAX y FTSE
dax_ftse_ret <- data.frame(diff(log(dax_ftse), lag=1))Ajustar un modelo de regresión no paramétrica NW entre los log-rendimientos de DAX en función de los log-rendimientos de FTSE. Obtener el ancho de banda óptimo utilizando un kernel gaussiano y utilizando un kernel uniforme. ¿Se puede decir que alguno de los modelos es mejor que otro en este caso?
# Gráfica de dispersión entre los log-rendimientos de DAX y FTSE
plot(dax_ftse_ret$FTSE,dax_ftse_ret$DAX)
# Ajustar un modelo de regresión no paramétrica NW entre los log-rendimientos de DAX en función de los log-rendimientos de FTSE.
# Ajuste del modelo NW GAUSSIANO
alpha <- seq(0.15, 0.99, by = 0.01)
a <- gcvplot(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern = "gauss", deg = 0, data = data8, alpha = alpha)
# Valor de alpha que minimiza los valores de a
alfa_star <- a$alpha[which.min(a$values)]
alfa_star## [1] 0.15modelo <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern = "gauss", deg = 0, data = data8, alpha = alfa_star)
summary(modelo)## Estimation type: Local Regression
##
## Call:
## locfit(formula = dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, data = data8,
## kern = "gauss", deg = 0, alpha = alfa_star)
##
## Number of data points: 1859
## Independent variables: dax_ftse_ret$FTSE
## Evaluation structure: Rectangular Tree
## Number of evaluation points: 36
## Degree of fit: 0
## Fitted Degrees of Freedom: 9.587b0 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="gauss", data = data8, deg = 0, alpha = 0.15)
b9 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="gauss", data = data8, deg = 0, alpha = 0.99)
m2 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="gauss", data = data8, deg = 0, alpha = alfa_star)
plot(dax_ftse_ret$DAX , col="darkgrey", main = "Nadaraya-Watson Gaussiano")
legend("bottomleft",
lty = 1,
col = c("red","blue","orange"),
legend=c(paste("bandwidth =", c( 0.15, 0.99, alfa_star))), bty = "n", lwd = 1.5)
lines(fitted(b0), col = "red", lwd = 1.5)
lines(fitted(b9), col = "blue", lwd = 1.5)
lines(fitted(m2), col = "orange", lwd = 3)
#--------------------------------------------
# Ajuste del modelo NW UNIFORME
alpha <- seq(0.15, 0.99, by = 0.01)
a <- gcvplot(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern = "rect", deg = 0, data = data8, alpha = alpha)
# Valor de alpha que minimiza los valores de a
alfa_star <- a$alpha[which.min(a$values)]
alfa_star## [1] 0.15modelo <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern = "rect", deg = 0, data = data8, alpha = alfa_star)
summary(modelo)## Estimation type: Local Regression
##
## Call:
## locfit(formula = dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, data = data8,
## kern = "rect", deg = 0, alpha = alfa_star)
##
## Number of data points: 1859
## Independent variables: dax_ftse_ret$FTSE
## Evaluation structure: Rectangular Tree
## Number of evaluation points: 36
## Degree of fit: 0
## Fitted Degrees of Freedom: 6.656b0 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="rect", data = data8, deg = 0, alpha = 0.15)
b9 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="rect", data = data8, deg = 0, alpha = 0.99)
m2 <- locfit(dax_ftse_ret$DAX ~ dax_ftse_ret$FTSE, kern="rect", data = data8, deg = 0, alpha = alfa_star)
plot(dax_ftse_ret$DAX , col="darkgrey", main = "Nadaraya-Watson Uniforme")
legend("bottomleft",
lty = 1,
col = c("red","blue","orange"),
legend=c(paste("bandwidth =", c( 0.15, 0.99, alfa_star))), bty = "n", lwd = 1.5)
lines(fitted(b0), col = "red", lwd = 1.5)
lines(fitted(b9), col = "blue", lwd = 1.5)
lines(fitted(m2), col = "orange", lwd = 3)
Obtener el ancho de banda óptimo utilizando un kernel gaussiano y
utilizando un kernel uniforme. ¿Se puede decir que alguno de los modelos
es mejor que otro en este caso?
No, porque en los dos se termina de la misma forma aunque se nota menor variabilidad en el de kernel uniforme